1 



Precision T Spectroscopy and Fundamental Parameters From NRQCD 

The NRQCD Collaboration 



On 



We present results from a high precision NRQCD simulation of the quenched T system at j3 = 6. We demon- 
strate a variety of important lattice techniques, including the perturbative improvement of actions, tadpole im- 
provement, and multicorrelated fits for extracting the spectrum of excited states. We present new determinations 
of a 3 (Mz) and M&, two fundamental parameters of the Standard Model. 
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1. Introduction 

In this paper we report new results from an 
accurate numerical simulation of the spectrum of 
the T family of mesons. (Results for the J/tp fam- 
ily are presented in (TJ.) Our simulation uses the 
NRQCD |Q action for the quarks, including all 
relativistic effects through 0(MbV 4 ) where Mt, is 
the 6-quark's mass and v its average velocity. We 
demonstrate a variety of important lattice tech- 
niques, including the perturbative improvement 
of actions, tadpole improvement, and multicorre- 
lated fits for extracting the spectrum of excited 
states. Furthermore, the high statistics and small 
systematic errors of our results allow us to ex- 
tract accurate values for two of the fundamen- 
tal parameters of the Standard Model: the mass 
of the 6-quark, and the strong coupling constant 

The 6-quarks in T's are quite nonrelativistic 
(v 2 is about 0.1). We exploit this in NRQCD 
by replacing the Dirac action for the quarks with 
a Shrodinger action. The Schrodinger theory 
is computationally much easier to solve because 
it can be treated as an initial-value problem, 
rather than a boundary- value problem. Relativis- 
tic effects are systematically introduced, order- 
by-order in v 2 , as corrections to the nonrelativis- 
tic action These terms are quite similar in 
form to the corrections that remove finite-lattice- 
spacing errors. We include both types of correc- 
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tion in our simulations. Consequently the dom- 
inant source of systematic error in our results is 
the gluon action — we use configurations produced 
with the standard Wilson action, and no light- 
quark vacuum polarization. (Finite volume errors 
are negligible since T's are much smaller than or- 
dinary mesons.) 

A complication in using improved actions such 
as ours is that each of the correction terms has its 
own coupling constant. These new coupling con- 
stants must be computed somehow. In principle 
they can be computed using weak-coupling per- 
turbation theory |^,^]. Our simulations demon- 
strate that this is also the case in practice, pro- 
vided tadpole-improved perturbation theory is 
employed. Indeed (tadpole-improved) tree-level 
perturbation theory seems quite sufficient for 
most aspects of T physics, at least at (3 = 6; and 
work has begun on the first-order corrections Q . 
This perturbative control over the effective action 
means that there are really only two parameters 
in the quark action, the mass and the charge, just 
as in the continuum theory. Thus our simula- 
tions are truly calculations from first principles, 
unlike calculations based upon a QCD-motivated 
phenomenological model like the quark potential 
model. 

In what follows, we first summarize the details 
of the simulation and fitting procedures. We then 
focus upon results for the spectrum, a -1 , Mb, and 
a-^(M z )- Finally we comment upon the impli- 
cations of our work concerning improved actions 
and simulations on coarse grids. 
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2. The Simulation 

The NRQCD quark lagrangian we used in our 
simulations was [B| 
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where n = 2, i?o is the nonrelativistic kinetic en- 
ergy operator, 
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and SH is the leading relativistic and finite- 
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Here A and A^ 2 ' are the simple gauge-covariant 
lattice derivative and laplacian, while A^ 4 ) is a 
lattice version of the continuum operator Y] Df. 
We used the standard cloverleaf operators for the 
chromo-electric and magnetic fields, E and B. 
The entire action was tadpole improved by divid- 
ing every link operator by uq = (^TrC/piaq) 1 / 4 . 
Potential models indicate that corrections be- 
yond SH contribute only of order 5-10 MeV to 
T energies. 

The only parameter in the NRQCD action is 
the bare mass of the quark, M" . We did a com- 
plete simulation for each of three masses: 2/a, 
1.8/a, and 1.71/a. Our analysis indicates that 
the last of these is the closest to the real mass 
(see below). However the statistical analysis for 
this mass is not yet complete and so most of the 
results we quote are for oMj = 1.8. Except where 
noted otherwise, the difference has negligible ef- 
fect. 

Note that meson energies in NRQCD are non- 
relativistic energies and do not include the rest 
mass energy. To obtain the rest mass of the T, 
for example, we measured the momentum depen- 
dence of its energy, fitting it to a form 

^(p)=^ R .(T) + ^£ m + ... (4) 



where the kinetic mass Mki n is the physical mass 
of the meson (see below). 

Our quark propagators were computed using 
an ensemble of quenched configurations obtained 
from the Staggered Collaboration ||. This en- 
semble contains 105 16 3 x 24 configurations, gen- 
erated with an unimproved Wilson action at J3 = 
6.0 and gauge fixed to Couloumb gauge. Potential 
models indicate that the errors due to quenching 
and to finite-et errors in the gluon action could 
be order 40 MeV in T energies, making these the 
most important systematic errors. 

We performed all of our measurements from 8 
different origins separated by L/2 in the initial 
time slice, and treated the propagators at these 
origins as statistically independent. At each value 
of the bare mass, our spectrum measurements re- 
quired a total of 17 heavy quark propagators, each 
with 3 spin-color components at the source and 6 
at the sink, using about 17 GFlops-hour of CPU 
time, i.e. about 800 inversions per GFlops-hour. 

When measuring meson propagators, we used 
gauge non-invariant smearing at both source and 
sink. With fast fourier transforms, all source- 
sink combinations can be measured at a compu- 
tational cost of one quark-propagator inversion 
per source smearing function Q. We used 15, 
2S, 35, IP, and 2P wavefunctions computed in 
the quark potential model (with the Richardson 
potential) as our smearing functions. In addition 
we used local sources for both 5 and P states. 
We also created smearing functions for 5-states 
with small, nonzero momenta. For each allowed 
set of quantum numbers, all source-sink smearing 
combinations were measured, resulting in a 4 x 4 
matrix of correlation functions for 5 states, 3x3 
for P states, and 2 x 2 for momentum states. 

3. Fitting 

We used a variety of fitting strategies to obtain 
the spectrum from our meson propagators. For 
the spectrum of spin-averaged 5 and P states we 
tried two sorts of multiexponential fit to multiple 
propagators. In the first we fit multiple expo- 
nentials to the set of propagators with smeared 
sources and a local sink: for example, the ground 
and first excited state of the T were determined 
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by simultaneously fitting three exponentials to 
the 1/ and 21 propagators, where 11 indicates 
IS smearing at the source and local smearing 
at the sink, and 21 indicates 2S smearing at the 
source and local smearing at the sink. We dis- 
carded results from the highest energy state in 
such fits, since this is the most susceptible to bi- 
asing due to higher states. These are the propa- 
gators for which our statistics are the best; how- 
ever, the statistical quality of the results degrades 
markedly for the excited states. 

The second approach we used for spin-averaged 
spectra was to fit the full matrix of propagators 
formed by taking every combination of nonlo- 
cal smearing function at the source and at the 
sink: for example, the T and its first two ex- 
cited states were determined by simultaneously 
fitting three exponentials to the 11, 12, 13, 21, 
22, 23, 31, 32, and 33 propagators formed us- 
ing quark-model wavefunctions for the first three 
S-states as smearing functions. The statistics 
were poorer for these smeared-smeared propaga- 
tors compared with the smeared-local propaga- 
tors discussed above; but the fitting is highly 
overdetermined, and consequently significantly 
more robust. This approach gave the best results 
for highly excited states; it also avoids spurious 
plateaus [Q. 

For splittings between two highly correlated 
states we found that the best method is to fit 
a single exponential to the jackknifed ratio of the 
two propagators for the two states. This approach 
is useful for computing spin splittings and the ki- 
netic mass of the T. 

Statistical uncertainties for the fitting param- 
eters were estimated by varying the parameters 
until Sx 2 = 1- In several cases we checked this 
estimate using a bootstrap analysis. The two esti- 
mates agreed in each case tested, suggesting that 
the statistical fluctuations in our propagators are 
at least approximately gaussian. 

Our fitting protocol follows closely the tech- 
niques discussed in ||, including the use of sin- 
gular value decomposition to allow the inversion 
of the variance matrix with finite precision arith- 
metic. In general, we had very few problems in 
obtaining successful fits, even for the second ex- 
cited T state. This is due to the high statistics 
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Figure 1 . NRQCD simulation results for the spec- 
trum of the T( 3 Si) and hbi}P\) and their radial 
excitations. Experimental values (dashed lines) 
are indicated for the S'-states, and for the spin- 
average of the P-states. The energy zero for the 
simulation results is adjusted to give the correct 
mass to the T. 



in our data set, and to our simultaneous fits us- 
ing many smearing functions. In particular, it 
was important to include a smearing function for 
every state we were attempting to fit. 

4. The Spectrum and a -1 

In Figure [l], we present our simulation results 
for the 3 Si and 1 P\ spectra, together with the 
corresponding experimental values from [[)). We 
use a -1 = 2.4 GeV and Monte Carlo data with 
aM® — 1.8. Because the hi, has not been ob- 
served, the experimental values quoted for the 
1 P\ states are actually the spin averages of the \b 
states; strong theoretical arguments, supported 
by our measurements, indicate that the two are 
the same. 

In Figure @, we present our results for the spin 
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Figure 2. Simulation results for the spin structure 
of the lowest lying P-state in the T family. The 
dashed lines are the experimental values for the 
triplet states, and the experimental spin average 
of all states for the singlet (hb). 



splittings of the P-wave ground state. Again 
a- 1 = 2.4 GeV and aM° = 1.8. The zero of 
energy is set to the spin average of the XbO-P)- 
The error bars are only for statistical errors. We 
expect systematic errors in the P-state energies 
of order 5 MeV, which is comparable to the sta- 
tistical errors shown. 

We determined the inverse lattice spacing by 
fitting the simulated spectrum to the experimen- 
tally determined spectrum. Rather than choos- 
ing a -1 to make, say, the Xb(lP) — T(15) split- 
ting correct, we used a bootstrap analysis to per- 
form a correlated fit for a -1 to the whole spec- 
trum. Specifically, treating the 840 origins in our 
ensemble as statistically independent, we gener- 
ated 20 bootstrap ensembles, each containing the 
meson propagators from 840 origins (with rep- 
etition). We then extracted the spectrum from 
each of these ensembles, giving us an ensemble 
of 20 spectra. Choosing various states for inclu- 
sion in the fit, we did a 1-parameter correlated fit 
of our spectra to the experimental spectrum to 
determine a -1 . This procedure was carried out 
twice, i.e. with two different sets of 20 bootstrap 
ensembles. 
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2.35(9) 
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Table 1 

Bootstrap results for a 



The results of our a -1 fits for different groups 
of states are summarized in Table |l|. The split- 
tings for 5-states are between 3 5i states, while 
spin averages are used for the P-states. The spin- 
splittings used are for the lowest lying P-state. 
For each set of states fit, we tabulate the values 
of a -1 and the associated goodness of fit Q for 
both sets of 20 spectra. (Q typically lies between 
0.1 and 0.9 for a good fit.) These fits indicate 
that 

a- 1 = 2.4(1) GeV. (5) 

Because of quenching, our a -1 should not agree 
with those obtained by matching simulation re- 
sults for observables that are more infrared, like 
the string tension, to experiment; we expect such 
determinations to yield smaller values of a" 1 . 

The success of our global fits of the T spectrum 
show that our simulation is accurately modeling 
the general features of the T physics. There are 
certainly systematic errors in our simulation re- 
sults, but these do not degrade the fits because 
they are generally smaller than the statistical er- 
rors. Only a small improvement in statistics is 
needed before these systematic effects will become 
apparent. Indeed there is already a slight indica- 
tion (2a) that the 25-15 and 1P-15 splittings 
are inconsistent, the first being too large and the 
second too small. This is precisely the effect ex- 
pected due to quenching of the gauge fields. Such 
discrepancies will be useful when comparable sim- 
ulations with unquenched configurations are be- 
gun. Similar effects are expected in the spin split- 
tings, particularly for 5-states. 
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aM» 


aS NR (T) 


aE Q 


M b 


6.0 


1.71 


0.455(1) 


0.21 


4.69 




1.80 


0.451(1) 


0.21 


4.70 




2.00 


0.444(1) 


0.22 


4.73 



Table 2 

Simulation and perturbative results used in the 
first method for determining Mb. Values for Mb 
are in GeV. 



P aM° aM kin (T) Z m I M kin (T) M b 
6.0 1.71 3.95(5) 1.15 9.5(4) 4.72 
1.80 4.10(5) 1.15 9.8(4) 4.76 
2.00 4.45(5) 1.14 10.7(4) 4.82 

Table 3 

Simulation and perturbative results used in the 
second method for determining M b . Values for 
M kin (T) and M b are in GeV. 



5. M b determination 

The &-quark's mass M b is an important funda- 
mental parameter of the Standard Model. Our 
simulation results lead to two independent deter- 
minations of this mass (ie, of the pole mass) []i~2[ . 
These are among the most accurate and reliable 
of all determinations to date. 

The first procedure for computing M b uses 
simulation results for the nonrelativistic en- 
ergy B N r(T) of the T (Eq. (|)). The quark mass 
is given by 

M b = \ (My - a- 1 (a£ NR (T) - 2aE )) , (6) 

where My = 9.46 GeV is the experimentally mea- 
sured mass of the T, and Eq is the nonrelativis- 
tic energy of a p = 6-quark in NRQCD. The 
quantity £'nr(T) — 2Eq can be thought of as the 
effective binding energy of the meson. The quark 
energy Eq is an ultraviolet divergent quantity and 
thus it can be computed using weak-coupling per- 
turbation theory 

aE = b a v {q ) (l + 0{a v )) (7) 

Here ay is the strong coupling constant as defined 
in (l(| , and qo ~ 1/a. In Table || we present the 
simulation results for -Bnr, and the correspond- 
ing Eq's from perturbation theory. From these 
results we conclude that the pole mass of the b- 
quark is M b = 4.7(1) GeV. The major source of 
uncertainty in this determination is from the two- 
loop corrections to E . These could be 10-20% of 
Eq, or 1-2% of M b . Uncertainties due to finite lat- 
tice spacing, quenching, tuning of the bare quark 
mass, statistics, and a -1 are all less than 1%. 
Our second procedure for determining M b is 



to tune the bare quark mass M 6 ° until the ki- 
netic mass of the T, as computed in the simu- 
lation (Eq. ([|)), agrees with the measured mass 
of the T. Then the pole mass of the quark is 
M b = Z m where renormalization constant Z m 
is computed using perturbation theory Q : 

Z m = l + b m a v (q m ) + 0(a v ) (8) 

where scale q rn ~ 1/a. To reduce the sensitivity 
of the result to the value of the lattice spacing and 
to the bare quark mass, we rewrite the expression 
for the pole mass as 

M b = Z m My "f) (9) 
aM kin (T) 

where My = 9.46 GeV is the T's experimental 
mass, while M^ in (T) is its mass as determined 
from our simulation with bare mass M®. Our 
results are summarized in Table |3|. These indi- 
cate that our best estimate for bare quark mass 
is = 1.7(l)/a, and therefore the pole mass is 
again M b — 4.7(1) GeV. Here the main source of 
uncertainty is again perturbative: the two-loop 
corrections to Z m could shift M b by 2-3%. Un- 
certainties due to a -1 , tuning the bare mass, etc. 
are of order 1% or less. 

Our two determinations of the 6-quark mass are 
very different, and yet they both yield a pole mass 
of M b = 4.7(1) GeV. The complete agreement 
between the two methods is a strong indication 
of the validity of each, and more generally of the 
lattice QCD techniques upon which they rely. 

6. a^(Mz) Determination 

The strong coupling constant is completely 
specified by two lattice quantities: the inverse 
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logWn 


0.152 


log W12 


0.152 


log W13 


0.152 


log W 2 2 


0.153 


log W 23 


0.154 


log W 33 


0.159 



Table 4 

The strong coupling constant as determined from 
Monte Carlo results for various Wilson loops. 



lattice spacing a -1 , and the coupling at a scale 
3.41/a. The inverse lattice spacing has been dis- 
cussed above. NRQCD simulations of the T sys- 
tem give a -1 = 2.4(1) GeV, where the errors are 
dominated by statistics. 

The coupling constant is obtained by fitting the 
Monte Carlo value of log(-|Tr[/ p i aq ) to the pertur- 
bative expansion fio|| 

— -ay(3.41/a){l-1.19ay}. (10) 

This procedure, in effect, defines ay. The re- 
liability of perturbation theory can be assessed 
by computing ay from a variety of Wilson loops, 
each with its own perturbative expansion (cor- 
rected to order ay). We have done this and the 
results are shown in Table ||. Loops with small ar- 
eas are less susceptible to nonperturbative effects 
and so we conclude that ay (3.41/a) = 0.152(1). 

Unfortunately we have rif — light-quark 
flavours in our simulation, while T physics is con- 
trolled by a theory with rif = 3. Thus a correc- 
tion is needed. The correction is illustrated by 
Figure [| which shows plots of the running cou- 
pling constants for the n/ = and rif = 3 the- 
ories that give the same T physics. An n/ = 
simulation gives rif = 3 results at scale g* if 

<$V) = «£V)- (11) 

For quarkonium the relevant g* is the mean 
momentum-transfer in the potential, which is 
computed from the expectation value of the po- 



tential between wavefunctions for the state of in- 
terest: 

Quark-model wavefunctions for Is, 2s, and 
Ip states imply g* is 0.75(25) GeV for T's. To 
correct for the wrong rif, we use the two-loop 
perturbative beta function to make the connec- 
tions: 

a { y } (3.41/a) 0^(3 GeV) 

I T (13) 

4°V) = 4 3 V) 

The perturbative beta function is probably re- 
liable down to scales of order 0.7-0.8 GeV (see 
Figure ||). We obtain a value for the unquenched 
coupling of 

a [ y ] (3 GeV) = 0.262(12) (14) 

where the uncertainty in q* is the dominant 
source of error. 

It is customary to quote values for the coupling 
constant in terms of a^(Mz)- The standard 
MS coupling is related to ay by 

a m (q) = a v (e 5 / 6 q) {l + 2a v /7r + 0(a 2 v )} . (15) 

Converting to MS and running the scale up to the 
mass of the Z meson, we obtain finally: 

= °- 112 ( 4 )' ( 16 ) 

which compares well with other lattice determina- 
tions plf , as well as with 0.110-0.125 as obtained 
from various high-energy experiments. The main 
sources of error in our result are: 



n / = 0^3 =^> 0.002 

ay — > =>■ 0.002 

(3.41/a) =S> 0.001 

a- 1 0.001 



All of these errors can be significantly reduced in 
the immediate future. 

7. Perturbative Improvement 

The lowest-order (or Schrodinger) NRQCD ac- 
tion accounts for the bulk of T and ip physics. 
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Figure 3. ^^(ijj for the n/ =0 and ra/ = 3 
theories having the same T physics. The two-loop 
perturbative beta function was used to generate 
the curves. The data points are values of oty) (q) 
determined from Monte Carlo measurements of 
various small Wilson loops and Creutz ratios of 
small loops. These show that two-loop perturba- 
tion theory is reliable over the entire range of q'a 
shown. 



However high precision results require an im- 
proved action. The key ingredients for an im- 
proved action are: 

• a tree-level perturbative calculation of the 
leading correction terms in the action 

(Eq. !)); 

• tadpole improvement of the correction 
terms, where every link is replaced by 

Up -» U^/uo, (17) 

and Uq is the Monte Carlo value of 
(^Trt/piaq); this step is essential for avoid- 



ing large renormalizations of the tree-level 
corrections; 

• a calculation of C(ay(l/a)) renormaliza- 
tions of the corrections (if very high pre- 
cision is needed). 

We have implemented the first two of these im- 
provements in our NRQCD simulations, and they 
have been strikingly successful, as is illustrated by 
the three examples listed in Table ||. 

The first example is the determination of the 
T mass from its dispersion relation (Eq. (^)). 
This is compared with the "exact" result obtained 
from the relation 

aMy = 2 (Z m oMf - aE Q ) + aE NR (T), (18) 

where renormalizations Z m and clEq are com- 
puted perturbatively, as discussed above (Eqs. (Q) 
and (||)). Both the addition of correction terms to 
the action and tadpole improvement are essential 
for accurate results. The most important correc- 
tion terms in the action for this quantity are the 
spin independent 0(v 2 7 a,a 2 ) terms. (Note that 
the kinetic mass of a composite particle equals 
the sum of the masses of the constituents in a 
nonrclativistic theory; the binding energy comes 
in only with relativistic corrections.) 

The second example is the spin-splitting be- 
tween the \b2 and Xbo P -state in the T family. 
This splitting vanishes in the lowest order the- 
ory since that theory is spin independent. Tad- 
pole improvement is essential here since the dom- 
inant operator that contributes to this splitting, 
%j)Hcr ■ D x g"Eil>, involves the cloverleaf operator 
for the chromoelectric field. Again only the fully 
corrected theory works well; it gives results in ex- 
cellent agreement with experiment. 

The third example is the spin-splitting between 
the ip and r\ c mesons jjj. The situation here 
is quite analogous to that for the P-state split- 
tings, except that here the dominant operator is 
t/^<t • g'Btp. Again correction terms and tadpole 
improvement are essential. The fully corrected 
simulation gives an excellent result; the 20% dis- 
crepancy between this result and experiment is 
consistent with the expected size of ay correc- 
tions, uncertainties in the mass, 0{v A ) operators, 
and effects due to quenching. 
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Quantity 


SH = 


"o = l 


Corrected 


Exact 




3.60(5) 


5.00(5) 


4.15(5) 


4.10(9) 


Xbl - XbO 





25(7) 


51(7) 


53 MeV 









94(1) 


117 MeV 



Table 5 

NRQCD predictions compared with exact results. Simulation results are presented for the theory without 
0(a,v) corrections (SH = 0), the theory without tadpole improvement (uq = 1), and for the fully 
corrected theory. 



These results are compelling evidence that per- 
turbative improvement of actions, when com- 
bined with tadpole-improved perturbation the- 
ory, is a highly effective procedure. 
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8. Conclusions 

We have presented early results from our con- 
tinuing program of high-precision analyses of 
heavy-quark mesons using NRQCD. We have 
shown that lattice simulations can accurately ac- 
count for the structure of the T spectrum up 
through the T(3S f ), and including spin structure. 
And we have used these results to make new de- 
terminations of the strong coupling constant and 
the 6-quark's mass. 

Our results are possibly the most thorough 
QCD tests to date of perturbative improvement 
for lattice actions. They underscore the reliability 
of perturbation theory, the utility of tree- level im- 
provement, and the critical importance of tadpole 
improvement (without it, the improvements tend 
to be much too small). The success of the per- 
turbative improvement for NRQCD strongly sug- 
gests that it will work for gluons and light quarks, 
particularly since the quarkonium states we ex- 
amine are 3-5 times smaller than light hadrons. 
(Note that 0{a 2 ) errors for the T at (5 = 6 are 
only of order a few percent, even though the me- 
son's radius is only two lattice spacings.) Sim- 
ulations of light hadrons with improved actions 
should permit precision work on lattices that are 
much coarser than those commonly used today. 
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